/*
 Copyright 2013--Present JMM_PROGNAME
 
 This file is distributed under the terms of the JMM_PROGNAME License.
 
 You should have received a copy of the JMM_PROGNAME License.
 If not, see <JMM_PROGNAME WEBSITE>.
*/
// CREATED    : 7/15/2015
// LAST UPDATE: 10/1/2015

#include "statsxx/optimization/bracket_methods/BracketMethod.hpp"

// STL
#include <cstdlib>    // std::abs()
#include <functional> // std::function<>, std::bind(), std::placeholders
#include <tuple>      // std::tuple<>, std::make_tuple()
#include <utility>    // std::swap()

// jScience
#include "jScience/linalg.hpp" // Vector<>, dot_product()
#include "jScience/physics/consts.hpp" // golden_ratio


inline BracketMethod::BracketMethod()
{
    this->max_stp = 100.0;
}

inline BracketMethod::~BracketMethod() {};


//========================================================================
//========================================================================
//
// NAME: std::tuple<double,
//                  double,
//                  double,
//                  double,
//                  double,
//                  double> BracketMethod::bracket(std::function<double(const double)> f,
//                                                 double x1,
//                                                 double x2)
//
// DESC:
//
// NOTES:
//    i. This subroutine ALWAYS returns a bracketed minimum.
//
//========================================================================
//========================================================================
inline std::tuple<
                  double,
                  double,
                  double,
                  double,
                  double,
                  double
                  > BracketMethod::bracket(
                                           std::function<double(const double)> f,
                                           double x1,
                                           double x2
                                           )
{
    double x3;
    
    double f1;
    double f2;
    double f3;
    
    f1 = f(x1);
    f2 = f(x2);
    
    if(f2 > f1)
    {
        std::swap(x1, x2);
        std::swap(f1, f2);
    }
    
    x3 = x2 + golden_ratio*(x2 - x1);
    f3 = f(x3);
    
    while(f3 < f2)
    {
        // use a parabolic fit between x1, x2, x3 to estimate the location of the minimum
        double t1 = (x2 - x1)*(f2 - f3);
        double t2 = (x2 - x3)*(f2 - f1);
        
        double u = x2 + ((x2 - x1)*t1 - (x2 - x3)*t2)/(2.0*(t2 - t1));
        double ulim = x2 + max_stp*(x3 - x2);
        
        double fu;
        
        // minimum between x2 and x3
        if( ((x2 - u)*(u - x3)) > 0.0 )
        {
            fu = f(u);
   
            // a minimum was found between x2 and x3
            if(fu < f3)
            {
                x1 = x2;
                x2 = u;
                
                f1 = f2;
                f2 = fu;
                
                break;
            }
            // ... minimum between x2 and u
            else if(fu > f2)
            {
                x3 = u;

                f3 = fu;
                
                break;
            }
            // ... parabolic fit was not useful
            else
            {
                // ... advance by the golden ratio
                u = x3 + golden_ratio*(x3 - x2);
                fu = f(u);
            }
        }
        // ... between x3 and ulim
        else if( ((x3 - u)*(u - ulim)) > 0.0 )
        {
            fu = f(u);

            // decreased, advance by the golden ratio
            if( fu < f3 )
            {
                x2 = x3;
                x3 = u;
                u = x3 + golden_ratio*(x3 - x2);
                
                f2 = f3;
                f3 = fu;
                fu = f(u);
            }
        }
        // ... beyond limit
        else if( ((u - ulim)*(ulim - x3)) >= 0.0 )
        {
            u = ulim;
        
            fu = f(u);
        
            // decreased, so advance by golden ratio
            if( fu < f3 )
            {
                x2 = x3;
                x3 = u;
                u = x3 + golden_ratio*(x3 - x2);
                
                f2 = f3;
                f3 = fu;
                fu = f(u);
            }
        }
        // ... reject parabolic fit, and advance by golden ratio
        else
        {
            u = x3 + golden_ratio*(x3 - x2);
            fu = f(u);
        }
        
        // shift the three points and continue
        x1 = x2;
        x2 = x3;
        x3 = u;
        
        f1 = f2;
        f2 = f3;
        f3 = fu;
    } // end while(f2 > f3)
    
    return std::make_tuple(x1, x2, x3, f1, f2, f3);
}